function [SpendRes,SpendIrr,DurRes,DurIrr,Prob,ExpSpend,x,F,rho,rhon,kappa,thetae,func,LHS,RHS,thetastar,thetahatind,star,starstar,nstar,nstarstar] = Sim(mu,sigma,xmin,xmax,N,beta,delta)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x        = linspace(xmin,xmax,N);
DEL      = (xmax-xmin)/(N-1);
ExpTheta = exp(mu+0.5*sigma^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1) Creating cdf and pdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

F     = cdf('lognormal',x,mu,sigma);
f     = pdf('lognormal',x,mu,sigma);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2) Calculate Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Q = 1 - F - x.*f*(1-beta);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3) Find argmin Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

thetahat = exp(mu+sigma^2/(1-beta));    %Analytical solution for lognormals
thetahatind = (thetahat-xmin)/DEL+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4) Define g and U
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computes flexible government spending and associated U(g)

g   = x./(x+beta*delta/(1-delta)*ExpTheta);
U   = log(g);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 5) Define rho
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computes theta^{**} (index j) given theta^* (index i)
%From equation 49

rho = zeros(1,floor(thetahatind));  %Works on indices, returns indices
j=N;

for i=1:length(rho)
    while sum(((Q(i:j-1)+Q(i+1:j))/2-Q(j))*DEL)<0
        j=j-1;
    end
    rho(i) = j;
end

if rho(end) == floor(thetahatind)
    rho(end) = rho(end)+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 6) Define rhon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computes theta_n^{**} (index j) given theta_n^* (index i)
%From equation 50

rhon = zeros(1,N);          %Works on indices, returns indices
j = 1;

for i=N:-1:ceil(thetahatind)
    while sum(((Q(j:i-1)+Q(j+1:i))/2-Q(j))*DEL)<0
        j=j+1;
    end
    rhon(i) = j;
end

if rhon(ceil(thetahatind)) == ceil(thetahatind)
    rhon(ceil(thetahatind)) = rhon(ceil(thetahatind))-1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 7) Define kappa
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computes theta_n^* (index j) given theta^* (index i)
%From equations 18 and 19

kappa = zeros(1,floor(thetahatind));  %Works on indices, returns indices
j=N;

for i=1:length(kappa)
    temp = sum(((U(i:rho(i)-1)+U(i+1:rho(i)))/2-U(i))*DEL);
    while j>thetahatind+1 && sum((U(j)-(U(rhon(j):j-1)+U(rhon(j)+1:j))/2)*DEL)>temp
        j=j-1;
    end
    kappa(i) = j;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 8) Find thetahat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computes, given theta^* (index i), the binding value of limited commitment
%constraint and the value of punishment. 
%Finds the theta^* that equalizes them
%From equation 20

func = zeros(1,floor(thetahatind));
thetastar = 0;
LHS = 0;
RHS = 0;
thetae = 0;
for i=1:floor(thetahatind)
    if rho(i)<N
        func(i) = sum(((U(i:rho(i)-1)+U(i+1:rho(i)))/2-U(i))*DEL) + delta/(1-delta)*(sum((Q(i:rho(i)-1)+Q(i+1:rho(i)))/2.*((U(i:rho(i)-1)+U(i+1:rho(i)))/2-U(i))*DEL)-sum((Q(rhon(kappa(i)):kappa(i)-1)+Q(rhon(kappa(i))+1:kappa(i)))/2.*((U(rhon(kappa(i)):kappa(i)-1)+U(rhon(kappa(i))+1:kappa(i)))/2-U(kappa(i)))*DEL));
        LHS(i) =  sum(((U(i:rho(i)-1)+U(i+1:rho(i)))/2-U(i))*DEL);
        RHS(i) = -delta/(1-delta)*(sum((Q(i:rho(i)-1)+Q(i+1:rho(i)))/2.*((U(i:rho(i)-1)+U(i+1:rho(i)))/2-U(i))*DEL)-sum((Q(rhon(kappa(i)):kappa(i)-1)+Q(rhon(kappa(i))+1:kappa(i)))/2.*((U(rhon(kappa(i)):kappa(i)-1)+U(rhon(kappa(i))+1:kappa(i)))/2-U(kappa(i)))*DEL));
    else
        thetae=i;
    end
end

thetastar = find(func<0,1);

if isempty(thetastar)
    thetastar = floor(thetahatind);
end

if thetastar == thetae+1
    thetastar = thetae;
end

%Expected spending from equations 12 and 14
SpendRes    = sum(g(1:thetastar-1).*f(1:thetastar-1)*DEL)+ sum(g(thetastar).*f(thetastar:rho(thetastar)-1)*DEL)+sum(g(rho(thetastar):end).*f(rho(thetastar):end)*DEL);
SpendIrr    = sum(g(1:rhon(kappa(thetastar))-1).*f(1:rhon(kappa(thetastar))-1)*DEL)+ sum(g(kappa(thetastar)).*f(rhon(kappa(thetastar)):kappa(thetastar)-1)*DEL)+sum(g(kappa(thetastar):end).*f(kappa(thetastar):end)*DEL);

%Expected duration of fiscal regimes
DurRes      = 1/(1-F(rho(thetastar)));
DurIrr      = 1/(1-F(rhon(kappa(thetastar))));

%Probability of fiscal regime, expected spending, variance
Prob        = (1-F(rhon(kappa(thetastar))))/(2-F(rhon(kappa(thetastar)))-F(rho(thetastar)));
ExpSpend    = Prob*SpendRes + (1-Prob)*SpendIrr;

%Equilibrium values of thetas
star = x(thetastar);
starstar = x(rho(thetastar));
nstar = x(kappa(thetastar));
nstarstar = x(rhon(kappa(thetastar)));
